* Mar 2015 Adam Sacarny - always run xtreg to generate ehat
* Jan 2012 Adam Sacarny
*! 1.0.2 5 Sep 2008 Austin Nichols
* 1.0.2  5 Sep 2008 fixes bug that ignores i() and absorb()
* 1.0.1 28 Feb 2008 makes abs option absorb and exits with error in no-regressor case
* 1.0.0 10 Feb 2008 Austin Nichols
program fese_adam, sortpreserve
	syntax varlist(numeric) [ if ] , ///
		HOmo(name) [HEtero(name) Ehat(varname)]
	
	* deal with grouping variable
	local group :char _dta[iis]
	
	* deal with varlist (LHS and RHS vars)
	if "`:word 2 of `varlist''"=="" {
		display as error "need at least 1 LHS and 1 RHS var"
		error 198
	}
	
	* deal with collinearity... bluntly
	* FIXME: this doesn't seem to deal with collinearity the only appears after
	* sweeping out the fixed effects...
	_rmdcoll `varlist'
	if (r(k_omitted) > 0) {
		display as error "collinearity in RHS variables"
		error 198
	}
		
	gettoken y x: varlist
	
	* deal with if()
	marksample touse
	
	* placeholders for standard errors
	gen `homo' = .
	if ("`hetero'"!="") {
		gen `hetero' = .
	}

	if ("`ehat'"=="") {
		tempvar ehat
		qui xtreg `y' `x' if `touse', fe
		predict `ehat' , e
	}
	
	sort `group'

	if ("`hetero'"=="") {
		mata: fese_o("`y'","`x'","`ehat'","`group'","`touse'","`homo'")
	}
	else {
		mata: fese_oh("`y'","`x'","`ehat'","`group'","`touse'","`homo' `hetero'")
	}
	
end

/*

prog fese, sortpreserve
 else {
  syntax [varlist] [if] [in] [, s(str) i(varname) Absorb(varname) Oonly * ]
  if "`i'`absorb'"!="" {
   if "`i'"!="" & "`absorb'"!="" {
    di as err "May not specify both " as txt "i" as err " and " as txt "absorb" as err " options."
    error 198
    }
   else loc gp "`i'`absorb'"
   }
  else loc gp:char _dta[iis]
  if "`:word 2 of `varlist''"=="" {
    di as err "You must specify at least one RHS (explanatory) variable for " as smcl "{help fese}"  as err "."
    di as err "A model without RHS (explanatory) variables can be estimated by " as smcl "{help mean}"
    di as err "with the " as smcl "{help mean:over()}" as err " option for het-robust SEs, or by hand for i.i.d. errors e.g.:"
    di as txt " areg y, a(" char(96) ":char _dta[iis]')"
    di as txt " predict fe, d"
    di as txt " replace fe=fe+_b[_cons]"
    di as txt " predict resid, r"
    di as txt " generate sqres=resid^2"
    di as txt " egen N=count(sqres), by(" char(96) ":char _dta[iis]')"
    di as txt " summarize sqres, meanonly"
    di as txt " generate se=sqrt(r(mean)*e(N)/e(df_r)/N)"
    di as err "See also " as smcl `"{browse "http://www.stata.com/statalist/archive/2008-02/msg00027.html":this Statalist post}"' _c
    di as err " for an illustration using " as smcl "{help xtreg}" as err "." _n 
    di as err "Clustered SEs would be zero (within the limits of machine precision) in this case"
    di as err "(for clustering on the panel [iis] variable as assumed elsewhere)." _n _n
    error 198
    }
  loc tv:char _dta[tis]
  qui g `s'se=.
  la var `s'se "OLS SE"
  if "`oonly'"=="" {
    qui g `s'hrse=.
    la var `s'hrse "Het-robust SE"
    qui g `s'crse=.
    la var `s'crse "Cluster-robust SE"
    loc sv "`s'se `s'hrse `s'crse"
  }
  else   loc sv "`s'se"
  marksample touse
  areg `varlist' if `touse', absorb(`gp') `options'
  gettoken y x: varlist
  qui {
   predict `s'b if `touse', d
   replace `s'b=`s'b+_b[_cons]
   la var `s'b "Estimated FE"
   tempvar e
   predict `e' if `touse', r
   sort `gp' `tv', stable
  } 
  if "`oonly'"=="" mata: fese("`y'","`x'","`e'","`gp'","`touse'","`sv'")
  else mata: fese_o("`y'","`x'","`e'","`gp'","`touse'","`sv'")
 }
end

*/

version 9.2
mata:
 void fese(string scalar depvar, string scalar x, string scalar r, string scalar G, string scalar tousename, string scalar s)
 {
  st_view(y, ., tokens(depvar), tousename)
  st_view(X, ., tokens(x), tousename)
  st_view(e, ., tokens(r), tousename)
  st_view(gp, ., tokens(G), tousename)
  st_view(sv, ., tokens(s), tousename)
  info = panelsetup(gp, 1)
  Ng=rows(info)
  N=rows(X)
  k=cols(X)+Ng
  sse=sum(e:^2)
  Xt=J(0,cols(X),0)
  for (i=1; i<=rows(info); i++) {
    Xi = panelsubmatrix(X, i, info)
    Xt=Xt\Xi:-mean(Xi)
  }
  XtXt=cross(Xt,Xt)
  _invsym(XtXt)
  OV=J(0,1,0)
  RV=J(0,1,0)
  CV=J(0,1,0)
  ov=J(0,1,0)
  rv=J(0,1,0)
  cv=J(0,1,0)
  for (i=1; i<=rows(info); i++) {
   ei = panelsubmatrix(e, i, info)
   Ti =(info[i,2]-info[i,1]+1)
   eiei=ei*ei'
   ei2=diag(ei:^2)
   di=J(info[i,1]-1,1,0)\J(Ti,1,1)\J(N-info[i,2],1,0)
   KDi=cross(Xt',cross(XtXt,cross(X,di)))
   KDii=KDi[info[i,1]..info[i,2],1]
   ci=sum(eiei)-2*cross(rowsum(eiei),KDii)
   ri=sum(ei2)-2*cross(rowsum(ei2),KDii)
   oi=Ti-2*cross(rowsum(I(Ti)),KDii)
   for (j=1; j<=rows(info); j++) {
    ej = panelsubmatrix(e, j, info)
    ejej=ej*ej'
    ej2=diag(ej:^2)
    KDji=KDi[info[j,1]..info[j,2],1]
    ci=ci+cross(KDji,cross(ejej,KDji))
    ri=ri+cross(KDji,cross(ej2,KDji))
    oi=oi+cross(KDji,cross(I(rows(ej)),KDji))
    }
   ooi=sqrt(sse/(N-k)*oi)/Ti
   rri=sqrt(N/(N-k)*ri)/Ti
   cci=sqrt((N-1)/(N-k)*(Ng/(Ng-1))*ci)/Ti
   ov=ov\ooi
   rv=rv\rri
   cv=cv\cci
   for (t=1; t<=Ti; t++) {
       OV=OV\ooi
       RV=RV\rri
       CV=CV\cci
       }
  }
  st_matrix("ov",ov)
  st_matrix("rv",rv)
  st_matrix("cv",cv)
  sv[.,.]=OV,RV,CV
}
void fese_o(string scalar depvar, string scalar x, string scalar r, string scalar G, string scalar tousename, string scalar s) {

	st_view(y, ., tokens(depvar), tousename)
	st_view(X, ., tokens(x), tousename)
	st_view(e, ., tokens(r), tousename)
	st_view(gp, ., tokens(G), tousename)
	st_view(sv, ., tokens(s), tousename)
	
	info = panelsetup(gp, 1)

	Ng=rows(info)
	N=rows(X)
	k=cols(X)+Ng

	sse=sum(e:^2)

	printf("sigma^2: ")
	printf(strofreal(sse/(N-k)))
	printf("\n")

	Xt=J(N,cols(X),0)
	
	for (i=1; i<=rows(info); i++) {
		Xi = panelsubmatrix(X, i, info)
		Xt[|info[i,1],1 \ info[i,2],cols(X) |] = Xi:-mean(Xi)
	}

	XtXt=cross(Xt,Xt)
	_invsym(XtXt)

	// ov=J(Ng,1,0)

	for (i=1; i<=rows(info); i++) {
	
		//printf("Working on object %f\n",i);
		//displayflush();
	
		//if ( mod(i,round(Ng/100))==0 ) {
		//	printf("Progress %f percent\n",100*i/Ng)
		//}
	
		Ti =(info[i,2]-info[i,1]+1)

		// printf("Starting...\n");
		// displayflush();
		

		Xdi = colsum(panelsubmatrix(X, i, info))'
		Xtdi = colsum(panelsubmatrix(Xt, i, info))'

		XtXtXdi = cross(XtXt,Xdi)

		dKDi = cross(Xtdi,XtXtXdi)
		dKKDi = cross(Xdi,XtXtXdi)
		
		oi=Ti-2*dKDi+dKKDi
		ooi=sqrt(sse/(N-k)*oi)/Ti
		
		//printf("Starting AN...\n");
		//displayflush();
		
		//di=J(info[i,1]-1,1,0)\J(Ti,1,1)\J(N-info[i,2],1,0)
		
		//KDi_an = cross(Xt',cross(XtXt,cross(X,di)))
		//dKDi_an = cross(di,KDi_an)
		//dKKDi_an = cross(KDi_an,KDi_an)
		
		//oi_an = Ti-2*dKDi_an+dKKDi_an
		//ooi_an=sqrt(sse/(N-k)*oi_an)/Ti
		
		//printf("Finished.\n");
		//displayflush();

		//ooi
		//ooi_an
		
		// ov[i,1]=ooi
		
		sv[| info[i,1],1 \ info[i,2],1 |] = J(Ti,1,ooi)
		
	}

	st_matrix("ov",ov)
}
 void fese_h(string scalar depvar, string scalar x, string scalar r, string scalar G, string scalar tousename, string scalar s)
 {
  st_view(y, ., tokens(depvar), tousename)
  st_view(X, ., tokens(x), tousename)
  st_view(e, ., tokens(r), tousename)
  st_view(gp, ., tokens(G), tousename)
  st_view(sv, ., tokens(s), tousename)
  info = panelsetup(gp, 1)
  Ng=rows(info)
  N=rows(X)
  k=cols(X)+Ng
  Xt=J(0,cols(X),0)
  for (i=1; i<=rows(info); i++) {
    Xi = panelsubmatrix(X, i, info)
    Xt=Xt\Xi:-mean(Xi)
  }
  XtXt=cross(Xt,Xt)
  _invsym(XtXt)
  RV=J(0,1,0)
  rv=J(0,1,0)
  for (i=1; i<=rows(info); i++) {
   ei = panelsubmatrix(e, i, info)
   Ti =(info[i,2]-info[i,1]+1)
   ei2=diag(ei:^2)
   di=J(info[i,1]-1,1,0)\J(Ti,1,1)\J(N-info[i,2],1,0)
   KDi=cross(Xt',cross(XtXt,cross(X,di)))
   KDii=KDi[info[i,1]..info[i,2],1]
   ri=sum(ei2)-2*cross(rowsum(ei2),KDii)
   for (j=1; j<=rows(info); j++) {
    ej = panelsubmatrix(e, j, info)
    ej2=diag(ej:^2)
    KDji=KDi[info[j,1]..info[j,2],1]
    ri=ri+cross(KDji,cross(ej2,KDji))
    }
   rri=sqrt(N/(N-k)*ri)/Ti
   rv=rv\rri
   for (t=1; t<=Ti; t++) {
       RV=RV\rri
       }
  }
  st_matrix("rv",rv)
  sv[.,.]=RV
}

void fese_oh(string scalar depvar, string scalar x, string scalar r, string scalar G, string scalar tousename, string scalar s) {
	st_view(y, ., tokens(depvar), tousename)
	st_view(X, ., tokens(x), tousename)
	st_view(e, ., tokens(r), tousename)
	st_view(gp, ., tokens(G), tousename)
	st_view(sv, ., tokens(s), tousename)
	
	info = panelsetup(gp, 1)
	
	Ng=rows(info)
	N=rows(X)
	k=cols(X)+Ng
	
	sse=sum(e:^2)
	
	Xt=J(N,cols(X),0)
	
	for (i=1; i<=rows(info); i++) {
		Xi = panelsubmatrix(X, i, info)
		Xt[|info[i,1],1 \ info[i,2],cols(X) |] = Xi:-mean(Xi)
	}
	
	XtXt=cross(Xt,Xt)
	_invsym(XtXt)

	OV=J(N,1,0)
	RV=J(N,1,0)
	
	ov=J(Ng,1,0)
	rv=J(Ng,1,0)

	for (i=1; i<=rows(info); i++) {
	
		if ( mod(i,round(Ng/1000))==0 ) {
			printf("Progress %f percent\n",100*i/Ng)
		}
		
		ei = panelsubmatrix(e, i, info)
		eiei=ei*ei'
		ei2=diag(ei:^2)

		Ti =(info[i,2]-info[i,1]+1)
		
		di=J(info[i,1]-1,1,0)\J(Ti,1,1)\J(N-info[i,2],1,0)
		KDi=cross(Xt',cross(XtXt,cross(X,di)))
		KDii=KDi[info[i,1]..info[i,2],1]
		
		ri=sum(ei2)-2*cross(rowsum(ei2),KDii)
		oi=Ti-2*cross(rowsum(I(Ti)),KDii)
		
		for (j=1; j<=rows(info); j++) {
			
			ej = panelsubmatrix(e, j, info)
			ej2=diag(ej:^2)
			
			KDji=KDi[info[j,1]..info[j,2],1]
			
			ri=ri+cross(KDji,cross(ej2,KDji))
			oi=oi+cross(KDji,cross(I(rows(ej)),KDji))

		}
		
		ooi=sqrt(sse/(N-k)*oi)/Ti
		rri=sqrt(N/(N-k)*ri)/Ti
		
		ov[i,1]=ooi
		rv[i,1]=rri
		
		OV[| info[i,1],1 \ info[i,2],1 |] = J(Ti,1,ooi)
		RV[| info[i,1],1 \ info[i,2],1 |] = J(Ti,1,rri)

	}
	
	st_matrix("ov",ov)
	st_matrix("rv",rv)
	sv[.,.]=OV,RV
}
end
exit

